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Abstract 


The present manuscript was written in 1994 and was not published. It addresses the form that 
the quantum-mechanical current density must take in mesoscopic treatments of semiconductor 
heterostructures, in which the electron dispersion relations are non-parabolic and position depen¬ 
dent, rendering the textbook expressions inapplicable. The approach is to derive the continuity 
equation for the specific model under consideration, using generalizations of Green’s identities to 
higher-order derivatives and to discrete models of different topological structure. A new addendum 
addresses two issues of more current interest: the use of irregular meshes in discrete formulations, 
and the identification of the Heisenberg velocity operator to evaluate current density. It is demon¬ 
strated that on discrete domains the velocity operator fails to satisfy a sensible continuity equation, 
and therefore cannot be identified with the current density. 


PACS numbers: 73.40.-c, 71.25.-s 


* frensley@utdallas.edu 


1 






ORIGINAL 1994 ABSTRACT 


In semiconductor heterostructures the electron dispersion relations are non-parabolic and 
vary with position. The non-parabolicity is described by effective Hamiltonians which either 
include higher-order derivatives or couple several basis states. As a result, the current density 
operator is not simply related to the gradient. By generalizing Green’s identity to higher- 
order derivatives and to difference relations, the appropriate form of this operator is derived 
for all of the commonly-used band structure representations. 


I. INTRODUCTION 

The wave mechanics of semiconductor heterostructures is complicated by the fact that the 
electron dispersion relation (or energy band structure) is generally non-parabolic at modest 
energies and necessarily varies with position. Under such circumstances, the form of the 
current density operator is no longer simply a symmetrized gradient. 

The form of the particle current density operator J is clearly constrained by the group- 
velocity theorem [l|, so that the expectation value of J on a state of definite wavevector k 
is 

{k\J\k) = Vg{k\k) = dE/dk{k\k)/h. (1) 

This equation, together with the band energies and eigenstates, in principle determines the 
form of J. It is, however, much more convenient to directly derive J from the Hamiltonian 
for a given problem. One does so by evaluating the time derivative of the probability density: 

d'il)*%l)/dt = (l/?h)['0*(if'0) — (if'0*)'0]. (2) 

Green’s identity [fV^g — = V ■ {fVg — gVf)], or a generalization thereof, is then 

invoked to write the right-hand side of ([2]) as the divergence of the current density. 

Heterostructures are most often described at a “mesoscopic” level where the microscopic 
(smaller than the atomic diameter) behavior of the wavefunction can be factored out. In 
order to realistically describe the non-parabolic dispersion relation, the resulting effective 
Hamiltonian must be more elaborate than a simple Laplace operator, and Green’s iden¬ 
tity must be correspondingly be generalized to derive the form of J. The commonly-used 
mesoscopic models can be classified as effective-mass or tight-binding approaches. 
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II. THE CURRENT DENSITY IN SPECIFIC REPRESENTATIONS 


For the purposes of the present discussion, we will assume that we possess a hermitian 
effective Hamiltonian which is valid throughout the structure, including any interval con¬ 
taining an abrupt heterointerface. The issues thus presumed to have been resolved have 
traditionally been posed as the definition of matching conditions for the mesoscopic wave- 
function. These matching conditions are frequently “derived” from the continuity of J Ml. 


B- 


but this condition is not sufficient to uniquely determine them [51. The Hamiltonian and the 
matching conditions are equivalent pieces of information, in the sense that either one may be 
derived from the other. In contrast, the continuity of J follows solely from the hermiticity 
of if je]. 

For the sake of simplicity, only the current density in one dimension will be considered. 
Extension to three-dimensional structures is in all cases straightforward, if somewhat tedious. 


A. Wannier-Slater Effective Mass Theory 

In the approach to effective-mass theory proposed by Wannier and expounded by Slater 
js], the microscopic wavefunction 'tp is expanded as a linear combination of localized Wannier 
functions, each of which is centered within a different unit cell. The expansion coefficients Tj 
can be regarded as values of a discrete lattice function which is interpolated between lattice 
points by a continuous function ^(z), for which an effective-mass Schroedinger equation is 
derived. If the electron dispersion relation can be expanded as 

N 

E{k) = ( 3 ) 


n=0 


then the effective Hamiltonian is 


N 


n=0 


( 4 ) 


dz^ 

To derive the appropriate J for this Hamiltonian, we require a generalized Green’s identity: 


Qn Qn Qn Qn £ Q 


n-1 


dz^ dz^ 


dz^ dz^ dz 




j=0 


Q3 j Qn j 1 ^Qng 


Qjg Qn j 1 ^9”/ 


dzi dz'^-i-^ dz^ dzi dz^-^-^ dz 


( 5 ) 
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(This identity is readily proven by expanding the derivative on the right-hand side; the 
summation then becomes a telescoping series.) The value of current density is thus 

N n—1 

n=l j=0 




Ar. 


-Ar, 


( 6 ) 


dz^ " dz"^ dz^ dz'^~^~^ " dz"^ 

Some formulations of the matching conditions produce terms in the Hamiltonian of the form 

+ d‘^^^^/dz^^^^B2n+i). (7) 

Their contribution to the current density may be readily derived from another identity: 




d2n+lg g2n+lf ^ g 
J dz^n+l +3 Q^ 2 n+l - ^ 


2n 






n-j 
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j=o 


dzj ’ 


leading to contributions to the current density of 


2n 




i=o 


d^B2n+l'4’* d^i’* ^B2n+li’ 


( 8 ) 


( 9 ) 


dz^ dz‘^^~^ dz^ dz‘^^~^ 

It appears (based upon an examination of some low-order cases) that any apparently her- 
mitian differential operator {q] can be manipulated into an expression containing only terms 
of the form (jl]) and ([7]). 


B. Luttinger-Kohn Effective Mass Theory 


The Luttinger-Kohn approach to effective-mass theory 


1CH12| more conveniently includes 


the effects of several bands. In this scheme the microscopic wavefunction is decomposed into 
the k = 0 Bloch functions and a set of slowly-varying envelope functions Xm{z), m being a 
band index. Differing numbers of bands may be included, with perhaps the most general 
scheme being that derived by Bastard [n|. By regrouping the various terms of Bastard’s 


Hamiltonian with respect to the derivatives, it can be written in the form 
implied over repeated indices) 


13| (summations 


f) f) 7 f) f) 

-jrAlmTT - Mmjr + 
oz oz 2 oz oz 


Xm, (10) 

where A and C are hermitian matrices and B a general matrix. H, H, and C are indexed 
by the band label m, and are z-dependent in a heterostructure. Using the above identities, 
the current density is readily shown to be 

dXn 


= 7 


ih 


X*iA^ 


Im ■ 


dz 


% 

-^AimXm + + B^)lmXr, 


Similar expressions have been obtained by Altarelli and by Burt 12 1 


( 11 ) 
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C. Tight-Binding Theories 


In the tight-binding approach 1^, the wavefunction is expanded in terms of a set of 
localized states m = 1,..., M in each atomic layer j 




( 12 ) 




The coefficients Cjm can be thonght of as forming a block-structured vector c with vector 
elements Cj = [cji,... ,CjM]'^■ The Hamiltonian then becomes a block-structured matrix 
of which the diagonal blocks are hermitian and describe interactions within a 

plane and the off-diagonal blocks are not necessarily hermitian and 

describe the coupling between planes. If only nearest-neighbor interactions are included, 
HTB) jg block-tridiagonal [^, 17 1 ^ 0 only for j = i — l,i,i + 1). 

Because the tight-binding representation is intrinsically discrete, we need to modify some¬ 
what our concepts of probability and current density. A total probability density pi is as¬ 
sociated with each atomic plane i, and is equal to cjcj. The current density represents the 
flux between adjacent planes; we will write the flux between planes i and z-f 1 as Jj+ 1/2 |l8 |. 
Applying ([2]) to the tight-binding Hamiltonian and assuming only nearest-neighbor 

interactions, we get 


d I 1 


ih . 




I «TTT(J^0J _ T TTI.J-DJ _ t 

This can be written as a discrete continuity equation, 

dpr'/9t = 

if is identified as 


{./'™’)i+i/2 = ^ (c.tiHggc, - c;H<™>c.+,) , 

Here we see the machinery of Green’s identity operating in a discrete space. 


(13) 


(14) 


(16) 


A variation of the tight-binding scheme is the “Wannier Orbital Model” 19|], which draws 
upon the discrete form of the Wannier-Slater theory. Interactions with remote neighbors are 
included to fit the E{k) dispersion, and typically M = 1, so that only one band is modeled. 
Inserting such a Hamiltonian into (E]) leads to many terms, which cannot be associated with 
a particular position. Thus, the notion of a local current density disappears, due to the 
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direct interactions between remote sites. Instead, we may define an antisymmetric current 
matrix with elements 


= (i/R) , 


( 16 ) 


where {J)ij is the current fiowing out of site i into site j. The nonlocal continuity equation 
is then 

= (17) 


III. SUMMARY 


The form of the current density has been derived for all of the usual models of het¬ 
erostructure electronic states. While one can usually find a way to obtain correct answers in 
a manual calculation without knowledge of the general expressions for J presented here, they 
provide a useful check on the results. The use of these expressions becomes more necessary 
if one seeks to develop the machinery for automatic computations (numerical or symbolic) 
applicable to a wide variety of heterostructures. 

The reader will have noticed that, in conformance to the practice in quantum-mechanics 
texts, expressions for the expectation values of J have been presented, not expressions for 
the operator itself. It is very difficult to represent J as an ordinary quantum-mechanical 
operator. Actually, J is much more naturally expressed as a superoperator which acts upon 


a density operator 


20|. All of the results presented here may be derived much more elegantly 


in a superoperator formalism, at the cost of an unfamiliar notation. 
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IV. 2015 ADDENDUM 


A. Irregular Triangular Meshes 


The finite-element technique has become a popular way of solving problems that are de¬ 
scribed by partial differential equations. One of the chief characteristics of this technique is 
the use of geometrically and topologically irregular meshes. Discussions of this approach al¬ 
ways cast it as merely an approximation to the continuum problem, but a more sophisticated 
approach will seek to determine the relations that a discrete model exactly satisfies. In par¬ 
ticular, we will ask what is the form of the continuity equation implied by the discretization 
of Schroediner’s equation on such an irregular mesh? 

Let us assume that Schroedinger’s equation has been discretized on a two-dimensional 
irregular triangular mesh using the “cotangent formula” for the Laplacian 21|. 






(18) 


where Uij is the sum of two cotangents of angles contained in the two triangles which share 
the i-j edge, and zero if the points i and j are more remotely located. Also, Uij = ooj^i and 
Ai is the area of the cell enclosing point i. If we use this Laplacian in the Schroedinger 
equation, and construct the continuity equation in the usual way, we find: 

^ ^ ^ {tp* - ^p*) - PiPj - pji)] , 

j 

j 


Thus, the finite-element case has the same structure as the Wannier-Orbital case, with a 
different current density associated with each pair of coupled mesh points. Observe that we 
have to take into account the variation of area associated with each mesh point, to express 
the continuity relation in terms of the total probability within each cell. 

One potential complication with this formulation is that the coupling elements Uij are not 
guaranteed to be non-negative, and negative values are known to occur if oblique triangles 
are present in the mesh. This will produce currents that flow in the “wrong” direction, but 


such effects also occur in solutions of the classical diffusion equation 


22 ]. 
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B. Relation Between Cnrrent Density and the Velocity Operator 


Another procedure which has been used to derive the current density is to identify it with 
the velocity operator, as dehned by the Heisenberg equation of motion: 


J^v = (i/h) [H, z] . (20) 

For a continuum model, this will produce the same results as ([6]) and ([9]) after symmetrization 
of V with respect to its adjoint. (That is, transformation to an anti-commutator superoper¬ 
ator.) 

For discrete formulations, this procedure produces a subtle discrepancy with respect to 
the discrete continuity equation. Consider the simple discretization of the effective-mass 
Hamiltonian and position operator on a uniform mesh of spacing a: 

~ 2fn*a^ ’ 


Then, 


= 


ih 

m* 


H+i,j 




2a 


( 21 ) 


Thus, the use of the Heisenberg equation leads to a velocity whose values are coincident 
with the mesh points and which is dehned by a centered diference. This is also what 
would be obtained by application of fl20|) followed by discretization according to textbook 
recommendations. 

Applying the adjoint symmetrization, we hnd the explicit form: 

ih 


{v)i = 




( 22 ) 


4m* a 

The tight-binding approach described above is equally applicable to this case, but it produces 
a current density value which is associated with the interval between meshpoints (as the 
i -|- 1/2 index implies), yielding the explicit expression: 


^ 2^ 1 (23) 

which will exactly satisfy the continuity equation ffTTj) . These differing results are related 
by: 

{'^)i = 2 (('^)i-l/2 + {■J)i+l/2) ■ (24) 
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In a steady-state problem in one dimension, all of the cnrrents discussed above will 
be equal, and equal to those at any other location. In transient situations, however, the 
central-difference will not exactly satisfy any simple continuity equation. Maintaining the 
central-diference assumption, the net inflow of {v) will be: 

= ^ {{J)i-3/2 + {J)i-l/2 - {J)i+l/2 - {J)i+3/2) , 

= 2 "^ + 4 (('^)*- 3/2 - {J)i+3/2) ■ ( 25 ) 


This is clearly not going to lead to any sensible continuity equation, as it invokes more 
remote current components. Consequently we must conclude that the Heisenberg velocity 
cannot be identihed with the current density in discrete domains. 

The expected rebuttal to the argument that I have just made is that the discrepancies 
will disappear as we let the mesh spacing approach zero. That cannot be taken for granted 
when the central-difference approximation to the gradient is employed. Because the central 
difference produces an anomalous value of zero for the maximum spatial Fourier component 
k = Tr/a, the convergence to the continuum is not uniform. This is the origin of the long¬ 
standing problem of anomalous states in envelope-function models of quantum states in 
heterostructures, and the remedy to this problem is the use of hrst-order differences 231. 
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